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ABSTRACT 



Context, The study of Type la supernovae (SNIa) has lead to greatly improved insights into many fields in astrophysics, e.g. cosmol- 
ogy, and also into the metal enrichment of the universe. Although a theoretical explanation of the origin of these events is still lacking, 
there is a general consensus that SNIa are caused by the thermonuclear explosions of carbon/oxygen white dwarfs with masses near 
the Chandrasekhar mass. 

Aims. We investigate the potential contribution to the supernova Type la rate from the population of merging double carbon-oxygen 
white dwarfs. We aim to develope a model that fits the observed SNIa progenitors as well as the observed close double white dwarf 
population. We differentiate between two scenarios for the common envelope (CE) evolution; the or-formalism based on the energy 
equation and the y-formalism that is based on the angular momentum equation. In one model we apply the tr-formalism always. In 
the second model the y-formalism is applied, unless the binary contains a compact object or the CE is triggered by a tidal instability 
for which the a-formalism is used. 

Methods. The binary population synthesis code SeBa was used to evolve binary systems from the zero-age main sequence to the 
formation of double white dwarfs and subsequent mergers. SeBa has been thoroughly updated since the last publication of the content 
of the code. 

Results. The limited sample of observed double white dwarfs is better represented by the simulated population using the y-formalism 
for the first CE phase than the a-formalism. For both CE formalisms, we find that although the morphology of the simulated delay 
time distribution matches that of the observations within the errors, the normalisation and time-integrated rate per stellar mass are a 
factor ~ 7 — 12 lower than observed. Furthermore, the characteristics of the simulated populations of merging double carbon-oxygen 
white dwarfs are discussed and put in the context of alternative SNIa models for merging double white dwarfs. 

Key words, stars: binaries: close, stars: evolution, stars:white dwarf, supemovae: general 



1. Introduction 

Type la supernovae (SNIa) are one of the most energetic explo- 
sive events known. They have been of great importance in many 
fields, most notably as a tool in observational cosmology. They 
have been used very success fully as standard c a ndles on cosmo- 
logica l distance scales (e.g. iRiess et al.l IT998I : iPerlmutter et alj 
1999), owing to th e special prop erty of great uniformity in the 
light curves (e.g. Phillips 1993). The SNIa also strongly af- 
fect the Galactic chemical evolution throug h the expulsion of 
iron (e.g. Ivan den Bergh & Tammannlll991 ). Despite their sig- 
nificance Type la supernovae are still poorly understood theoret- 
ically. 

Supernovae Type la are generally thought to be caused by 
thermonuclear explosions of carbon/oxygen (CO) white dwarfs 
(WD s) with masses near the Chandrasekhar mass M c h ~ 1 AM G 
(e.g. lNomotoi ri982). Various progenitor scenarios have been 
proposed. The standard scenarios can be di vided into two 
schoo ls of thoughts: the single-degenerate (SD) (Wh elan & Ibenl 
19731) and double-deg enerate (DD) scenario (Webbink 1984; 
Iben & Tutukovll 19841) . In the SD scenario, a CO WD explodes 
as an SNIa if its mass approaches M c h through accretion from 



a non-degenerate companion. In the DD scenario, two CO WDs 
can produce an SN la while merging if their combined mass is 
larger than M^. 

However, observationally as well as theoretically, the ex- 
act nature of the SNIa progenitors remains unclear. The explo- 
sion mechanism is complex due to the interaction of hydrody- 
namics and nuclear reactions. Several models exist that vary 
for example between a detonation or deflagration disruption or 
vary between explosions at th e Chandrasekhar mass or sub- 
Chandrasekhar masses (see e.g. Hille brandt & Niemeverfl2000l 
for a review). It also remains unclear whether the DD and SD 
scenario both contribute to the SNIa rate or if one of the scenar- 
ios dominates. Both scenarios have problems in matching the- 
ories with observations. A serious concern about the DD sce- 
nario is whether the collapse of the remnant would lead to a su- 
per nova or to a neutron star through accretion-induced collapse 
("seelNomoto & Ibenl 1 9851: ISaio & Nomotoll985blPiersanti et al ' 
2003; lYoon et al.l l2007t iPakmor et all I201QL l2012t IShen et al 
2012 :). Although in the SD channel the models for the explo- 
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sion process need to be fine-tuned to reproduce the observed 
spectra and light curves, an SNIa like event is more easily re- 
produced in the simulations of the explosion process. One prob- 
lem with the SD scenario is that the white dwarfs should go 



1 



S. Toonen et al.: Supernova Type la progenitors from merging double white dwarfs: 



through a long phase of supersoft X-ray emission, although it 
is unclear if there are enough of these sources to account for 
the SNIa rate (see iDi Stefanoll2010t iGilfanov & Bogdaril20Tot 
lHachisu et al.l 120101) . Moreover archival data of known SNIa 
have not shown this em ission unambiguously, but there is ma 
be one case (see Vbss & Nelemansl 120081; iRoelofs et al.l 1200 
INielsen et aUboioF Furthermore, SNIa that take place more 
than a few 10 9 years after the starburst (see e.g. Maoz et al. 2010) 
are hard to create in this cha nnel (e.g.fYungelso n & Liviol l2000; 
IHan & Podsiadlowskill2004 . 

To use SNIa as proper standard candles, we need to know 
what SNIa are, when they happen and what their progenitors 
are. Therefore, we study the binary evolution of low- and inter- 
mediate mass stars. In a forthcoming paper (Bours, Toonen & 
Nelemans, in preparation) we study the SD-scenario by look- 
ing into the poorly understood physics of accretion onto white 
dwarfs. In this paper we focus on the DD scenario and the ef- 
fect of the as yet very uncertain phases of common envelope 
(CE) evolution on the double white dwarf (DWD) population. 
These DWD systems are interesting sources for studying var- 
ious phases of stellar evolution, in our case the CE evolution. 
Gravitational wave emission is also important as this affects 
the binary system by decre asing the orbital period and even- 
tually leading to a merger dKraft et all 1 1 962t1Petersl 1 1 964 . or 
possible a SNIa. The DWDs are expected to be t he dominant 
source dEvans et al.l 119871: iNelemans et al.l 1200 lab of gravita- 
tional waves for the fut ure space-born gravitational wave obser- 
vatories such as eLISA dAmaro-Seoane et"aT1l2012a| b|). 

We study the population of merging DWDs that might lead 
to a SNIa from a theoretical point of view. We incorporated re- 
sults from observations where possible. We use the population 
synthesis code SeBa for simulating the stellar and binary evo- 
lution of stellar systems that leads to close DWDs. In Sect. [2] 
we describe the code and the updates since the last publication 
of SeBa. A major influence on the merging double-degenerate 
population is the poorly understood C E phase dPaczynskill 1 976t 
Webbink 1984; Ne lemans et ai1l2000l) . We adopt two different 
models for the CE. In Sect. [3] we describe these models and their 
implications for the observations of close DWDs. In Sect. [4] we 
discuss the binary paths leading to SNIa for each model. The 
SNIa rates and time-integrated numbers are derived in Sect. [5] 
The properties of the population of merging DWDs are discussed 
in the context of the classical and alternative sub- and super- 
Chandrasekhar SNIa explosion models in Sect. [6] A discussion 
and conclusion follows in Sect. [7] 

2. SeBa - a fast stellar and binary evolution code 

We present an update to the software package SeBa 

dPortegies Zwart & Verbuni 119961; INelemans et aD l2001bl) for 
fast stellar and binary evolution computations. Stars are evolved 
from the zero-age main sequence (ZAMS) until remnant for- 
mation and beyond. Stars are parametrised by mass, ra- 
dius, luminosity, core mass, etc. as functions of time and 
initial mass. Mass loss from winds, which is substantial 
e.g. for massive stars and post main-sequence stars, is in- 
cluded. Binary interactions such as mass loss, mass trans- 
fer, angular momentum loss, CE, magnetic braking, and grav- 
itational radiation are take n into account with appropriate 
recipes at every timestep dPortegies Zwart & Verbuntl [l996; 
IPortegies Zwart & Yungelsonlll998l) . Following mass transfer in 
a binary, the donor may turn into a helium-burning star without 
hydrogen envelope. When the mass transfer leads to a merger 
between the binary stars, we estimate the resulting stellar prod- 



uct and follow the evolution further. Note that we do not solve 
the equations of stellar structure. The stellar tracks instead as- 
sume stellar models in hydrostatic equilibrium. When this is not 
the case, however the gas envelope surrounding the core may 
puff outward (see Appendix IA.2I for details on the formalism). 
In our simulation the mass transfer rate is calculated from the 
relevant timescales (see Appendix IA.31 > and not from than the 
stellar radii. Therefore binary evolution is not critically depen- 
dent on out-of-equilibrium parameter values. 

The philosophy of SeBa is to not a priori define the binary's 
evolution, but rather to determine this at runtime depending on 
the parameters of the stellar system. When more sophisticated 
models become available of processes that influence stellar evo- 
lution, these can be included, and the effect can be studied with- 
out altering the formalism of binary interactions. An example is 
the accretion efficiency onto the accretor star during mass trans- 
fer. Instead of prescribing a specific constant percentage of the 
transferred matter to be accreted (and the rest to be lost from the 
system), the efficiency depends on the properties of the accreting 
star, such as the thermal timescale, the radius and the Roche lobe 
of the accretor (see Appendix lA.2l for details). Another example 
is the stability of mass transfer. In our simulations the stability 
and rate of mass transfer are dependent on the reaction to mass 
change of the stellar radii and the corresponding Roche lobes. 
The advantage of this is that the (de)stabilisin g effect of non- 
conse rvativeness of stable mass transfer (see Soberma n et al.l 
1 1997b is taken into account automatically. There is no need to 
make the assumption in the stability calculation that stable mass 
trans fer is conservative, as with methods that depend on the mass 
ratio dHiellming & Webbinldl987tlTout et all 1 9971: iHurlev et al.l 
120021) . 

Since the last publication of the code content, many changes 
have been made. We briefly discuss the most important changes 
below, and provide more detail in Appendix [A] First, the wind 
mass loss prescriptions th at we implemented a re mostly based on 
the recommendations by iHurlev et al.l (120001) . The specific pre- 
scriptions for different types of stars are described in Appendix 
IA.1I Second, a summary of the treatment of accretion onto dif- 
ferent stars can be found in Appendix IA.2I The accretion pro- 
cedure previously used in SeBa is complemented with a pro- 
cedure for accretion from a hydrogen-poor star. We assume 
that for ordinary stars, helium-rich matter is accreted directly 
to the core of the star. The mass accretion process onto white 
dwarfs is updated with new efficiencies of mass retention on 
the surface of the white dwarf. For hydrogen acc retion we have 
the op tion t o choose between t h e effic iencies of lHachisu et al.l 
(2008f) and iPrialnik & Kovetzl d 19951). Helium re tention can 
be modelle d according to iKato & Hachisu (1999) (with up- 
dates from lHachisu et alJ 119991) or llben & Tutukovl d 19961). I n 
this research we used t he ef ficiencies of Hachisu et all (|2008), 
IKato & Hachisul d 19991) . and lHachisu et alU 1999). For a study 
of different retention efficiencies and the effect on the Supernova 
Type la rate using the new version of SeBa, see Bours, Toonen 
& Nelemans, in preparation. Third, the stability of mass trans- 
fer is based on the adiabatic and thermal response of the donor 
star to mass loss and the response of the Roche lobe. The adjust- 
ment of the Roche lobe is dependent on the mass transfer rate, 
which in turn sets the efficiency of accretion onto the accretor 
star, see Appendix lA.3l Fourth, regarding the stellar tracks, pre- 
viously, stellar evolution has been based on evolutionary tracks 
described by analytic formulae given bvlEggleton et al.l (l989, 
hereafter EFT) with updates from [Tout et all (1 19971) and he- 
lium s tar evoluti on as described b y Porteg ies Zwart & Verbuntl 
(Il996l) based on llben & Tutukovl dl985l) . In the new version. 
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Fig. 1. Simulated population of visible double white dwarfs as a function of orbital period and mass of the brighter white dwarf. 
Left: the stellar evolution tracks according to EFT are used; right: HPT (using model ya, see Sect. [3]). The intensity of the grey scale 
corresponds to the density of objects on a linear scale. The sa me grey scale is use d for both plots. Obse rved binary white dwarf s 
are overplotted with fill e d circles. Thick points taken are f romlMarsh et al] (1201 ll). thinner points from Tovmassi an et al] j2 004); 
iNapiwotzki etal] d2005h : rKulkarni & van Kerkwiikl ( l20Toh : lBrown et al] d2010ll201 ll) ; lMarsh et al] d201 ll) : lKilic et al] d201 laldlbl) . 
see Sect. l2.1l for a discussion. 



the evolution of ordinar y stars and hydrogen-poor stars is based 
on iHurlev et all d2000l hereafter HPT). We do not adopt the 
HPT tracks for remnants, instead we maintain our prescrip - 
tion dPortegies Zwart & Verbunll 19961: iNelemans etaLll2001bh . 
which includes processes such as natal kick velocities to com- 
pact objects in supernovae explosions. 



2.1. Impact on the population of double white dwarfs 

Fig. Q] shows the visible close DWD population simulated by 
SeBa. On the left a simulation is shown of the previous ver- 
sion of SeBa that a.o. uses the EFT tracks and on the right we 
show the current version using the HPT tracks. Initial param- 
eters are distributed according to the distributions described in 
TableQ] Primary masses are drawn from 0.96-1 IM Q to include 
all stars that evolve into a white dwarf in a Hubble time. For 
the mass ratio and eccentricity we cover the full range 0-1, and 
the orbital separation out to IO 6 R . We assumed solar metallic- 
ity, unless specified otherwise. In the normalisation of the sim- 
ulation we assumed that primary masses lie in the range 0.1- 
1OOM . Our m ethod to estimate the v isible population of DWDs 
is described in Nele mans et al.l d2004h. in which the Galacti c star 
formation history is base d on iBoissier & Prantzosl d 1999b and 
WD cooling according to iHansenl dl999t) . We assume a mag- 
nitude limit of 21. In Fig.Q] the observed DWDs are overplotted 
with filled circles. The systems are d escribed by iMarshl d201 ll) 
and references therein, as we ll as iTovmassian et al.l ([2004 ) 
and iRodrfguez-Gil et al.l d2010l). Additionally, we included 1 9 
newly discov e red DWDs fro m lKulkarni & van Kerkwiikl ( 20 1 Oh ; 
irownj 
( 2011 



Table 1. Distributions of the initial binary parameter mass, mass 
ratio, orbital separation and eccentricity. 



y discov e red DWDs fro m lKulkarni &"va n Kerkwiik (2010); 
met al] (120101 120111): Marsh et all (120111) : lKilicetal.1 
laid) and iKilic et al.l (1201 lbk These new systems are dis- 



played with smaller circles and thinner lines to separate them 
from the previously found systems. We did this because the ob- 
servational biases are very different. The previously found sys- 
tems were selected from a magnitude-limited sample down to 
16-17 magnitude. The new systems are much fainter at about 20 
magnitude. Moreover, m ost of the new syst ems are discovered as 
part of the ELM survey (Bro wn et alj2010l) . This survey focuses 
on finding extremely low-mass white dwarfs from follow-up ob- 



Parameter 



Mass of single stars 
Mass of binary primaries 
Mass ratio 
Orbital separation 
Eccentricity 



Distribution 
Kroupa IMF"* 
Kroupa IMF (1) 
Flat distribution 

N{d)da x ar l da (2) 
Thermal distribution (3) 



Notes. ( " lKroupaetaII ( fl993h : (2) r AblfT983l) : <3) |Heggie ( 1975) 



servations of spectroscopically selected objects from the Sloan 
Digital Sky Survey. Therefore, the set of new systems is biased 
to lower masses. One should take this bias into account while 
comparing with the simulations and not take the combined set of 
obser ved systems as a rep resentative sample of the DWD popu- 
lation. |KiHc_eTay (|20lla]) showed in their Fig. 12 a visualisation 
of the population of visible DWDs simulated by SeBa, where 
this selection effect has been taken into account. 

The locations of the observed DWDs in Fig. Q] correspond 
reasonably well to the predictions of both models. The overall 
structure of the simulated populations from both models are sim- 
ilar. At masses of ~ O.5M and periods of 1-10 hr, there is a very 
pronounced region in the plot from the EFT tracks that seems to 
be missing in the HPT plot. However, this is not really the case. 
These systems mainly consist of one helium (He) WD and one 
CO WD. Masses of CO WDs span a wider range of values in 
the HPT tracks, which distributes the pronounced region in EFT 
over a larger region in mass and period in HPT. 

For a single burst of star formation the number of created 
DWDs within 13.5 Gyr and with an orbital period P < lOOOhr 
for the HPT and EFT stellar tracks is very similar; 6.9 ■ 10~ 3 
per M of created stars for both models. The time-integrated 
merger rate is 2.4 ■ 1O" 3 M ' for HPT and 3.2 ■ 1O" 3 M ' for 
EFT. The current merger rate in the Milky Way according to 
the HPT and EFT stellar tracks is very similar; 1.4 ■ 10~ 2 yr~' 
for HPT and 1 .2 • 10~ 2 yr~' for EFT, for which we have assumed 
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a star formation his t ory as in Nelema ns et al.l (|2004) based on 
iBoissier & Prantzosl(ll999l) . 

Classically, the population of double He dwarfs is thought to 
dominate in number over the other types of close DWDs. Using 
the EFT tracks for a single burst of star formation, we predict a 
percentage of [He-He, He-CO, CO-CO] = [60%, 17%, 21%] and 
a negligible number of DWDs containing oxygen/neon (ONe) 
dwarfs. For the HPT tracks, the percentage of double He dwarfs 
decreases to 38%. The population consists of [He-He, He-CO, 
CO-CO]=[38%, 27%, 33%] and 2% CO - ONe dwarfs. The de- 
crease in number of He WDs is caused by a difference in the 
stellar tracks related to helium i gnitio n under degenerate con- 
ditions. As shown bv lHan et al.l (|2002), degenerate stars do not 
ignite helium at a fixed core mass, but instead the core mass at 
helium ignition is a decreasing function of the ZAMS mass of 
the star. Taking this into account, more WDs in close binaries 
are labelled CO WDs. 



3. Two models for common envelope evolution 

Close DWDs are believed to encounter at least two phases of 
mass transfer in which one of the stars loses its hydrogen enve- 
lope. In at least one of these phases mass transfer from the evolv- 
ing more massi ve star to the less massive compa nion is dynam- 
ically unstable dPacz vnskill 1 976t IWebbink 1984) which leads to 
a common envelope. The core of the donor and companion spi- 
ral inward through the envelope, expelling the gaseous envelope 
around them. Because of the loss of significant amounts of mass 
and angular momentum, the CE phase plays an essential role in 
binary star evolution in particular the formation of short-period 
systems that contain a compact object. 

Despite of the importance of the CE phase and the enormous 
efforts of the community, all effort so far have not been success- 
ful in understanding the phenomenon. Several prescri ptions for 
CE ev olution have been proposed. The c-formalism (Webbink 
1 19841) is based on the conservation of orbital energy. The a- 
parameter describes the efficiency with which orbital energy is 
consumed to unbind the CE according to 



a(E, 



orb,init 



J orb,final 



(1) 



where E or \, is the orbital energy and E gr is the binding energy 
between the envelope mass M em and the mass of the donor M. 
E gr is often approximated by 



p — 



GMM e , 



AR 



(2) 



where R is the radius of the donor star and A depends on the 
structu re of the donor. We assume aA - 2. iNelemans et al.l 
(2000) deduced this value from reconstructing the last phase of 
mass transfer for 10 known DWDs using the unique core-mass 
- radius relation for giants. 

To explain the observed distribution of DWDs, 
INelemans et al.l (120001) proposed an alternative formalism. 
According to this y-formalism, mass transfer is unstable and 
non-conservative. The mass-loss reduces the angular momentum 
of the system in a linear way according to 



-Anil - -/final AM 

■ = y- 



Jir 



M + TO ' 



(3) 



where J mlt resp. 7fl na i is the angular momentum of the pre- and 
post-mass transfer binary respectively, and to is the mass of th e 
companion. We assumed y = 1.75, see INelemans et all 120 lb). 



We adopt two evolutionary models that differ in their treat- 
ment of the CE phase. In model aa the (^-formalism is used 
to determine the outcome of every CE. For model ya the 
y-prescription is applied unless the binary contains a com- 
pact object or the CE is triggered by a tidal instability (rather 
than dynamically unstable Roche lobe overflow, see Appendix 
IA.3I >. Typically, the second CE (with a giant donor and white 
dwarf companion) is described by the c-formalism, which 
gives consistent result s when compared with the observations 
(INelemans et al.ll2000h . If the first phase of mass transfer is un- 
stable, it typically evolves through a y-CE. In model ya and aa, 
if both stars are evolved when the CE develops, we assumed that 
both cores spiral-in (see Nelemans et al. 2001b). The envelopes 
are expelled according to 



J gr,*don 



J gr,*comp 



a(E, 



orb,init 



J orb, final 



(4) 



analogous to Eq[T] where Zs^+don represents the binding en- 
ergy of the envelope of the donor star and £g ri * CO mp of the com- 
panion star. 

The motivation for the alternative formalism is the large 
amount of angular momentum available in binaries with similar- 
mass objects. The physical mechanism behind the y-description 
remains unclear however. Interesting to note here is that recently 

Woods et a n (120101 120T1 suggested a new evolutionary model 

to create DWDs. It differs from standard assumptions in the first 
phase of mass transfer. These authors find that mass transfer be- 
tween a red giant and a main-sequence star can be stable and 
nonconservative. The effect on the orbit is a modest widening, 
with a result alike to the y-description. 

3.1. Impact on the population of double white dwarfs and 
type la supernova progenitors 

Fig. [2] shows the mass ratio of the visible population (see Sect. 
I2.ll ) of DWDs versus the orbital period according to model ya 
and aa. Overlayed with filled circles are the observed popula- 
tions. For systems for which only a lower limit to the mass of the 
companion is known, we show a plausible range of mass ratios 
of that system with an arrow. The arrow is drawn starting from 
the maximum mass ratio, which corresponds to an inclination 
of 90 degrees. It extends to a companion mass that corresponds 
to an inclination of 41 degrees. Within this range of inclinations 
there is a 75% probability that the actual mass ratio lies along 
the arrow. The filled circles overplotted on the arrow indicate 
the mass ratio for the median for random orientations, i.e., 60 
degrees. 

Using model aa, the DWDs cluster around a mass ratio of 
q ~ 0.5, while model ya shows a wider range in mass ratio. This 
agrees better with the observed binaries. The different mass ratio 
distributions are inherent to the models and only slightly depen- 
dent on the CE efficiency. This is because in the first CE phase, 
the y-CE allows for widening or very mild shrinkage of the or- 
bit, whereas in the a-prescription the orbit will always shrink. 
The resulting orbital separation determines when the secondary 
will fill its Roche lobe, and the corresponding core mass of the 
secondary, which determines the mass ratio distribution of the 
prospective DWD. 

In Fig.[3]the population of observed and simulated DWDs are 
shown as a function of combined mass of the two WDs for the 
two models of CE evolution. The left upper corner bounded by 
the dotted and dashed lines contains SNIa progenitors. In Fig. 
[3] there are two systems that have a probability to fall in this 
region. These systems are the planetary nebulae nuclei with WD 
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log P(hr) 

(a) 



log P(hr) 

(b) 



Fig. 2. Simulated population of visible double white dwarfs as a function of orbital period and mass ratio, where mass ratio is defined 
as the mass of the brighter white dwarf divided by that of the dimmer white dwarf. In the left model ya is used, in the right model 
aa. The intensity of the grey scale corresponds to the density of objects on a linear scale. The same grey scale is used for both plots. 
Observed binary white dwarfs are overplotted with filled circles, see Fig.Q]for references. 



companions TS 01 (PN G135.9+55.91(lTovmassianet alJl2004t 
lNapiwotzkietalJl2005l;[Tovmassianetal.ll2010l) and V458 Vul 
(Rod riguez-Gil et al]l2010l) . An immediate precursor of a DWD 
that is possibly a progenitor candidate for a SN la via the DD 
channel has also been observed ; a subdwarf with a white dwarf 
comp anion, KPD 1930+2752 dMaxted et all 120001: iGeier et alJ 
120071) . 

In our model of the visible population of DWDs (see Sect. 
12. U , the percentage of merging DWDs with a total mass ex- 
ceeding the Chandrasekhar is 1.2% for model ya and 4.3% for 
model aa. Including only double CO WDs, the percentage is 0.9 
and 2.9%, respectively. Because the number of observed close 
DWDs until today is low, we do not expect to o bserve many 
SNIa progenitors (see also iNelemans et al]|2001bl) . Therefore a 
comparison of the SNIa progenitors with population synthesis 
by a statistical approach is unfortunately not yet possible. We 
find it important to compare the observed close DWD population 
with the simulated one, since these systems go through similar 
evolutions and are strongly influenced by the same processes. 
Although the observed population mostly consists of He DWDs 
and He - CO DWDs instead of CO DWDs required for SNIa 
progenitors, at this time the population of all close DWDs are 
the closest related systems that are visible in bulk. 

4. Evolutionary paths to supernova type la from the 
double degenerate channel 

In this section we discuss the most common binary scenarios that 
leads to a potential supernova type la in the DD channel. We as- 
sume that every merger of two carbon/oxygen white dwarf with 
a mass exceeding 1 .4M will lead to a supernova. The contribu- 
tion of merging systems that contains a helium white dwarf that 
surpasses the Chandrasekhar mass is negligible. In the canonical 
scenario a DWD is formed through two consecutive CEs. This 
we label the 'common envelope channel'. In accordance with 
iMennekens et alJ ([2010), we find that there are other channels 
that can lead to a SNIa as well. We find that the common en- 
velope scenario can account for less than half of the supernova 
progenitors in a single burst of star formation, 34% for model ya 
and 45% for model aa. We distinguish between three scenarios 
labelled 'common envelope', 'stable mass transfer' and 'forma- 



tion reversal' . The names of the first two tracks refer to the first 
phase of mass transfer, whereas 'formation reversal' applies to 
the reversed order in which the two white dwarf are formed, see 
Sect. 14.3b . The stable mass transfer channel accounts for 52% 
and 32% assuming model ya and aa, respectively, for a single 
burst of star formation. The formation reversal channel accounts 
for a lower percentage of all SNIa, 14% for model ya and 23% 
for model aa for a single starburst. Note that the importance 
of the stable mass transfer channel strongly depends on the as- 
sumed amount of mass loss and angular momentum loss. 

In population synthesis studies all known information 
about binary evolution is combined, and different evolution- 
ary paths emerge out of these quite naturally. As noted by 
Menneke ns et alJ (12010). the significant contribution to the SNIa 
rate from other channels than the common envelope channel 
complicates the use of analytical formalisms for determining the 
distribution of SNIa delay times. The SNIa delay time of a binary 
is the time of the SNIa since the formation of the system. This 
is commonly used to compare observational and synthetic rates 
to constrain different physical scenarios (e. g. Yungelson & Liviol 
l2000tlRuiter et al.ll2 009: Men nekens et al.ll2010i see also Sect. 151 
in this paper). During a CE phase the companion is assumed to 
be hardly affected e.g. by accretion. If this is not the case, as in 
stable mass transfer, the assumption that the formation timescale 
of the DWD is approximately the main-sequence lifetime of the 
least massive component is not valid any more. Furthermore, 
the in-spiral timescale from DWD formation to merger due to 
gravitational waves is strongly dependent on the orbital sepa- 
ration at DWD formation. This can be very different for sys- 
tems that undergo stable mass transfer instead of a CE evolu- 
tion. Concluding, the delay time, which is the sum of the DWD 
formation and in-spiral timescale can be significantly different 
when these tracks are not properly taken into account. 



4. 1 . Common envelope channel 

In the canonical path, both stars lose their hydrogen envelopes 
through two consecutive common envelopes. An example of a 
typical evolution is shown in Fig. [4] In this example two zero-age 
main-sequence stars of 6M and 4M are in an orbit of 125 days. 
When the initially more massive star (hereafter primary) ascends 
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Fig. 3. Simulated population of visible double white dwarfs as a function of orbital period and the combined mass of the two dwarfs. 
On the left the common envelope is parametrised according to model ya, on the right according to model aa (see Sect. [3}. The 
intensity of the grey scale corresponds to the density of objects on a linear scale. The same grey scale is used for both plots. Observed 
binary white dwarfs are overplotted with filled circles, see Fig.Q]for references. The Chandrasekhar mass limit is indicated with the 
dotted line. The dashed line roughly demarks the region in which systems merge within a Hubble time. Systems located to the left 
of the dashed line and above the dotted line are supernova type la progenitors in the standard picture. 
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4.01 30.03 
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4.01 31.17 CE - WD + G 
0.55 0.17 WD + WD 



Fig. 4. Evolutionary track for the merger of two car- 
bon/oxygen white dwarfs of a combined mass that exceeds the 
Chandrasekhar mass. In this scenario the first phase of mass 
transfer is dynamically unstable which results in a common en- 
velope. In this figure we show a representative example of model 
ya, see Sect. [3] ZAMS stands for zero age main sequence and 
CE for common envelope. G is a giant star, MS a main-sequence 
star, He MS a helium MS and WD a white dwarf. 



pens, mass transfer is usually dynamically stable, but the effect 
on the orbit is small. 

A variation of this evolution can occur when the secondary 
has reached the giant stages of its evolution when the primary 
fills its Roche lobe. This happens for systems of nearly equal 
masses. We assume both stars lose their envelope in the CE 
phase according to Eq. |4] in which the orbit is severely de- 
creased. This variation contributes 23% of the systems in the 
common envelope channel for model ya and 10% for model aa. 

With the a-CE prescription it is likely to have another vari- 
ation on the evolution, in which the primary becomes a white 
dwarf immediately after the first phase of mass transfer. This 
can happen when the primary fills its Roche lobe very late on 
the asymptotic giant branch when the star experiences thermal 
pulses. These systems have initial periods that are a factor 5 lger 
than in the standard CE channel using model aa. This subchan- 
nel contributes 43% to the CE channel for model aa. When us- 
ing the ya-model for the CE, the contribution from this subchan- 
nel is 20%. However, these systems are not formed through a 
standard y-CE because the orbit does not shrink severely enough 
to obtain a significant contribution. Instead these systems are 
formed through a double-CE as described by Eq|4] For model 
aa the double-CE mechanism is important in only 18% of the 
43%. 



the giant branch, it fills its Roche lobe and a CE commences. The 
primary loses its hydrogen envelope, but does not become a WD 
immediately and a helium star is born. The primary becomes a 
white dwarf of about solar mass. When the initially less massive 
star (hereafter secondary) evolves off the main sequence and its 
radius significantly increases, another common envelope occurs. 
As a result, the orbit shrinks. The secondary evolves further as 
a helium star without a hydrogen envelope until it eventually 
turns into a white dwarf. For model aa the orbit decreases more 
severely in the first phase of mass transfer. Therefore the initial 
periods in this channel are higher, by a factor ~ 1 .5 - 3 and the 
primaries are typically more evolved giants when the donors fill 
their Roche lobes. In this evolution channel both the primary and 
secondary can fill their Roche lobes as helium giants. If this hap- 



4.2. Stable mass transfer channel 

In this channel the initial masses of the stars and the initial orbits 
are smaller than for the common envelope channel. Typical val- 
ues are a primary mass of 5M G , a secondary mass of 3M Q and 
an orbital separation of 40R Q (assuming a circular orbit). The 
primary fills the Roche lobe as a Hertzsprung gap star and mass 
transfer occurs stably. Which fraction of transferred mass is ac- 
tually accreted by the secondary and how much is lost from the 
system depends on the mass and radius of the secondary and the 
secondary's Roche lobe (see Appendix lA.2l for more details). In 
Fig. |5h an example of a typical evolution is shown. When the 
secondary fills its Roche lobe, a CE commences. In this channel 
the tidal instability (see Appendix lA.3l ) is important. In one third 
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Fig. 5. Two evolutionary tracks for the merger of two carbon/oxygen white dwarfs of a combined mass exceeding the Chandrasekhar 
mass. In these scenarios the first phase of mass transfer is dynamically stable. On the left an example of the stable mass transfer 
channel is shown, on the right the formation reversal channel. In the latter scenario the first phase of mass transfer is dynamically 
stable which results in a low-mass helium-star with a long lifetime. The initially less massive star becomes the first formed white 
dwarf. The top and bottom parts of the figure have different scales due to a common envelope phase, denoted as CE in the figure. 
Abbreviations are as in Fig. [4] Additionally sRLOF stands for stable Roche lobe overflow, HG is a Hertzsprung-gap star and He G 
a helium giant. 



of the systems the CE occurs because of a tidal instability, the 
other part is caused by a dynamical instability. The secondary 
turns into a hydrogen-deficient helium-burning star in a system 
in which the period has decreased by one or two orders of mag- 
nitude. As in the previous channel, the primary and secondary 
can fill their Roche lobe as helium giants. If the primary fills 
its Roche lobe, mass transfer is usually dynamically stable and 
has little effect on the orbit. If the secondary fills its Roche lobe, 
mass transfer can be stable or unstable. In the example of Fig. [5^ 
when the secondary fills its Roche lobe again as it ascends the 
helium giant branch, the mass transfer is unstable and the orbital 
separation decreases by a factor ~ 5. 

4.3. Formation reversal channel 

We present a scenario in which in the first mass transfer a helium 
star (sdB star) is formed that becomes a white dwarf only after 
the companion has become a white dwarf.Q A typical example of 
an evolution like this is shown in Fig.|5j5. The first phase of mass 
transfer is stable, like the stable mass transfer track. However, 
the resulting helium stars in this channel have low masses in 
the range of 0.5-0. 8M and long lifetimes of ~ 10 8 yr. The first 
mass transfer occurs approximately conservatively. As a conse- 

1 Th is track is a close analogy of the track proposed by Sipio r et alJ 
(2004) regarding recycled pulsars. In the scenario proposed by these 
authors the end states of the two components are reversed, resulting in 
a neutron star that forms prior to a black hole. However, in our scenario 
the name 'formation reversal' applies to the evolutionary timescales of 
the primary and secondary. Although the primary first evolves off the 
main sequence, the secondary becomes a remnant first. 



quence, the subsequent evolution of the high-mass secondary 
(5-8M ) accelerates. When the secondary fills its Roche lobe, 
mass transfer is tidally unstable (see Appendix IA.3I >. The sec- 
ondary loses its hydrogen and helium envelope in two consecu- 
tive CEs and becomes a white dwarf. Subsequently, the original 
primary evolves of the helium main-sequence and becomes a 
white dwarf. 

To our knowledge, this track has not been studied in detail 
before. Therefore we evaluated this track by performing detailed 
numerical calculations using the ev b inary stellar-evolution 
code originally develop ed by Eggleton (Egglet on! 1 19711 1 19721 
I Yakut & E ggleton 2005, and r efere nces therein) and upda ted as 
described in lPols et alJ d 19951) and iGlebbeek et alJ 12008). The 
code solves the equations of stellar structure and evolution for 
the two components of a binary simultaneously. The simulation 
showed that indeed the evolution of the secondary can be accel- 
erated through accretion so that the secondary can stop helium 
burning prior to the primary. 

5. Delay time distribution 

One way to constrain the population of SNIa progenitors is 
through the delay time distribution (DTD), where the delay times 
is the time between the formation of the binary system and the 
SNIa event. In a simulation of a single burst of star formation 
the DTD gives the SNIa rate as a function of time after the star- 
burst. The DTD is linked to the nuclear timescales of the pro- 
genitors and the binary evolution timescales up to the merger. 
We assumed a 50% binary fraction and initial parameters are 
distributed according to the distributions described in Table [TJ 
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Fig. 6. Merger rate of double carbon/oxygen white dwarfs with 
a total mass above the Chandrasekhar mass as a function of 
delay time. Rates are in yr~' per 1O 1O M formed stellar mass 
of the parent galaxy. Solar metallicity (Z = 0.02) is assumed. 
Delay times are shown for two different prescriptions of the 
CE phase. In black we plot model ya and in grey model 
aa, see Sect. [3] Overplotted with black circles are the ob- 
served values of the SNIa rate of | Totani et al. (2008), Maoz et alJ 
d2010h. 



., .Ma oz & Badenes] d2010l) and iMaoz et all d2011l) (see 
Maoz & Mannuccill20 1 1 , for a review). For comparison the grey 



circles show the observations scaled down by a factor 10. 
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Fig. 7. Merger rate of double carbon/oxygen white dwarfs with a 
total mass above the Chandrasekhar mass as a function of delay 
time for a metallicity of 0.001. Rates are in yr~' per 1O 1O M per 
formed stellar mass of the parent galaxy. Delay times are shown 
for two different prescriptions of the CE phase. In black model 
we plot ya and in grey model aa, see Sect. [3] 



In Fig. [6] we compare the delay time distribution for the two 
different models of CE evolution. The sharp cut-off near 13.5 
Gyr is artificial, because evolution was only allowed to proceed 
for 13.5 Gyr. The delay time distribution shows that these merg- 
ers are expected to take place in young as well as old popula- 
tions. The peak in the supernova la rate is at ~150 Myr for both 
models. The median delay time is 0.7 Gyr for model aa and 1.0 
Gyr for model ya. The normalisations of the delay time distri- 



bution of model aa and ya are comparable. The time-integrated 
number of SNe la per unit formed stellar mass is 2.0 • 1O~ 4 M ' 
and 3.3 • lO^M" 1 for model ya and aa, respectively. From the 
Lick Observatory Supernova Search. lMaoz et all (1201 ll) inferred 
a value of 2.3 ±0.6 • 10~ 3 M Q ' , which is a factor 7-12 higher than 
the predictions from our models. 

The morphologies of the DTD of model ya and aa resemble 
each other in that they show a strong decline with delay time, 
although with a slightly different slope. The aa model shows 
higher rates at short delay times, whereas the rate for model 
ya shows higher rates at long delay times. This is because the 
a-CE causes a stronger decrease of the orbital separation than 
the y -CE in the fi r st pha s e of mass transfer. The observed rates 
fromlTotani et all d2008l). IMaoz et all d2010l), IMaoz & Badenesl 
d2010l) . and lMaoz et alJ ( 120 1 lh fsee lMaoz & Mannuccil201 ll for 
a review), shown in Fig. |6l are much higher than the predicted 
rates from both models. To compare the morphological shapes 
of the DTDs more easily, we scaled down the observations by a 
factor 10 in Fig.|6]in light grey. The shape of the observed DTD 
fits the synthetic DTDs well. At long delay times > 6 Gyr, the 
flattening of the DTD is better reproduced by the ycr-model. 

We have a last remark ab out Fig. [6] about the dat- 
apoint from IMaoz et al .] ( 120 Id) at 185Myr and a rate of 
0.165 yr~' (10 10 M o ) _1 . If this datapoint is true, it could indi- 
cate a steep rise of the delay time distribution at the shortest 
delay times. Neither model ya, or aa reproduces the steep rise 
indicated by this point. At short delay times the contribution to 
the SNIa rate from other channels might be significant, for in- 
stance the contrib u tion from he l ium d onor s in the SD channel . 
IWang et al.1 d2009h . iRuiter et alJ d2009h and lClaevs etall d2011h 
showed that at delay times of ~ lOOMyr, the DTD from he- 
lium donors in the SD channel peaks, although rates at this de- 
lay time vary between 10~ 4 - 10~ 2 yr~' (10 10 M Q ) . Hydrogen 
donors in the SD channel are a possible contributor to the SNIa 
rate as well, but there is a strong disagreement oyer the DTD 
from this channel, (see for example Nelemans et al. 2012, for an 
overview). In that paper it was shown that the simulated peaks of 
the DTDs lie anywhere between 0. 1 to 3 Gyr and the peak rates 
vary between 10~ 6 - 10 _2 yr _1 (lO lo Af ) -1 . 

If we do not assume an instantaneous burst of star formation, 
but instead convolve the DTD with a star formation rate, we can 
estimate the SNIa rate from double degenerates for spiral galax- 
ies like th e Milky Way. If we as sume a Ga lactic star formation 
rate as in Nelemans et al .1 d2004 based on Bois sier & Pra ntzos 
d!999h . model aa gives 8.3 • 10" 4 SNIa yr _1 . Model ya gives 
a Galactic rate of 5.8 • 10" 4 SNIa yr . The reason for the rela- 
tively high Galactic rate for model ya in comparison with model 
aa relative to the integrated rates is that the peak of star for- 
mation occurs at long delay times where the DTD of model 
ya dominates over model aa. The empirical SN la rate fr om 
Sbc-type galaxies like our own (Cappellaro &Turattd l2001) is 
4 + 2 • 10~ 3 yr -1 , which is a factor ~ 5 - 7 higher than the simu- 
lated rates. 

When convolving the DTD with a star formation history, one 
should also take into account a metallicity dependence of the 
stars. To study the effect of metallicity on the SNIa rate, we sim- 
ulated a delay time distribution from a single burst of stars of 
metallicity Z = 0.001. The important part of the DTDs in this 
respect are the long delay times because this is where the frac- 
tion of metal poor stars is highest. The DTD of model ya for 
Z = 0.001 is lower at long delay times and higher at short de- 
lay times than the same CE model for solar metallicities. The 
time-integrated number of SNe la per unit formed stellar mass 
is 2.8 ■ 1O~ 4 M ' for model ya. This is an increase of ~60% 
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with respect to Z = 0.02. The integrated rate for model aa for 
Z = 0.001 is 4.2 • 10~ 4 M-\ which is an increase of -30% with 
respect to Z = 0.02. The DTD for Z = 0.001 is roughly simi- 
lar in morphological shape to that for Z = 0.02, see Fig. [7] The 
DTDs of both metallicities for model aa are similar at long de- 
lay times, and consequently the effect on the Galactic SNIa rate 
is expected to be marginal. 

6. Population of merging double white dwarfs 

In this section we discuss the properties of the population of 
SNIa progenitors from merging DWD and place it in the con- 
text of recent studies of the SNIa explosion itself. Fig. [8] shows 
the combined mass of the system as a function of delay time 
for merging CO DWDs. It shows that for classical SNIa pro- 
genitors, the number of merging events decreases with time and 
that the number decreases faster with time for model aa than for 
model ya, as discussed in Sect. [5] Moreover, the figure shows 
that mergers near the Chandrasekhar mass are most common, 
inde pendent of delay time. 

iFrver et ail (1201 Ol) showed that if super-Chandrasekhar 
mergers of CO DWDs of ~ 2M produce thermonuclear ex- 
plosions, the light curves are broader than the observed SNIa 
sample. These authors argued that these mergers cannot dom- 
inate the current SNIa sample. We find indeed in both models 
that mergers with combined masses ~ 2M are much less com- 
mon than mergers in systems with a combined mass near the 
Chandrase khar mass limit. 

Where IFrver et al.1 d2010h studied a mer g er of a 1 .2M Q CO 
WD with a O.9M CO WD, iPakmor et all (1201 Oh focused on 
mergers of nearly equal mass WDs. In their scenario both WDs 
are distorted in the merger process and the internal st r ucture 
of the merger remnant is quite different. Pak mor et al.l ([2010) 
argued that these mergers become hot enough to ignite car- 
bon burning if the WD masses exceed M > O.9M . They 
found that t hese systems re semble subluminous SNIa such as 
SN 1991bg. lLi~etaLl J2001b found 1991bg-like supern ovae ac- 
count for 16 ±6% of all SNIa. From an improved sample lEi et al.l 
(1201 lh found a percentage of 15.2j?. If we assume that 1991bg- 
like events account for 15% of all SNIa and the time-integ rated 
of all SNIa types is 2.3 + 0.6 ■ lO^M" 1 dMaoz et al.ll201ll) . the 
time-integrated rate should be 3.5 ± 0.9 • 10~ 4 M ' . When we in- 
clude the error of 6% on the fraction of 1991bg-like events with 
respect to all SNIa, the rate is 3.5 + 1.6- 1O~ 4 M '. If we assume in 
a simplistic way that all CO DWD mergers of q = M 2 /M\ > 0.92 
and Mi > 0.9 (where M\ is the most massive WD and M2 the 
least massive WD) would lead to a 1991bg-like event, the time- 
integrated number of events is 2.3 ■ 1O~ 5 M ' according to model 
ya and 1.8 ■ IO^Mq 1 assuming model aa. While the SNIa rate 
from the classical progenitors from model ya is comparable to 
that of model aa, the population of DWDs is very different. In 
Sect. 13. II we showed that the type of CE parametrisation intro- 
duces a bias in the mass ratio distribution of observed DWDs, 
which mostly consist of (double) He DWDs and He-CO DWDs. 
In Fig. [9] we show that this is also the case for the population of 
merging CO DWDs. Although the mass ratio distribution is not 
important for the sta ndard DP scenario, i t is important for the 
scenario proposed by IPakmor et al.l ( 2010[). If the s tandar d sce- 
nario and the scenario proposed by Pakmor et alJ (|2010) hold, 
SN 1991bg-like events are more common in model ya. We have 
to make a side remark on the expected delay times of this sce- 
nario. The median delay times are 180 Myr for model ya and 
150 Myr for model aa. The timescales are short because gen- 
erally, more massive WDs have more massive progenitor stars, 



whose evolutionary timescales are short compared to those of 
less massive stars. Observations show, however, that sublumi- 
nous SNIa are associated with old stellar populations ~ 5-12Gyr 
(lHowellll2001h . 

In our simulations the majority of merging CO DWDs have 
combined masses below the Chandrasekhar mass, see Fig. [8] 
Sub-Chandrasekhar models have long been propose d in order 
to rais e total number of SNIa to match observations. ISim et al.l 
(2010) found that if sub-Chandrasekhar WDs can be detonated, 
especially in the range ~ 1.0-1 .2M , the explosions match sev- 
eral observed properties of SNIa reasonably well. In the hypo- 
thetical situation that all double CO WDs that merge lead to an 
SNIa event, the integrated rate is 8.3 ■ 1O~ 4 M ' for model ya and 
9.3- IO^Mq 1 for model aa. This is still a factor 3 lower than the 
observed rate of 2.3 + 0.6 • 10 3 M ' dMaoz et al.ll20TTh . Only if 
we assume that all mergers between a CO WD and a CO or He 
WD lead to an SNIa, the rates of 1.6 ■ IQ^M^ and 2.1 ■ lO^Af" 1 
for model ya and aa, respectively, match the observed rate. 

The challenge for sub-Chandrasekhar models is how to deto- 
na te the white dwarf. A scen ario for this was recently suggested 
by van KerkwijketaT](l2010l) . In this scenario two CO WDs with 
nearly equal masses merge. The merger remnant itself is too cold 
and insufficiently dense to produce an SNIa by its elf, as noted by 
IPakmor et al.l d2010l) . Ivan Kerkwiik et all (|2010) proposed that 
accretion of the thick disk that surrounds the remnant leads to 
an SNIa through compressional heating. If we simplistically as- 
sume that every merger of a double CO DWD with q > 0.8 and 
M\ + M2 < 1 .4 leads to an SNIa, the time-integrated number per 
unit formed mass is l^-lO^M" 1 for model ya and 8.8- 10 5 M Q l 
for model aa. Relaxing the condition of Mi + M2 < 1 .4 to all 
masses, 2.3 ■ IO^Mq 1 for model ya and 1.9-10~ 4 M-' for model 
aa. As in the scenario proposed by Pakmor et af] d2010l) . when 
a scenario is biased to merging systems of high-mass ratio, the 
relative contribution from this scenario in the ya model is higher 
than the aa model. 



7. Conclusion and discussion 

We studied the population of SNIa progenitors from merg- 
ing double CO WDs with a combined mass exceeding the 
Chandrasekhar mass, the so-called DD progenitors. We consid- 
ered two prescriptions of the CE phase. The CE evolution is 
a crucial ingredient in the formation of close double degener- 
ate compact objects, but the process itself is still poorly under- 
stood. The first model assumes the a-formalism for all CE. The 
second model is a combination of the a-formalism and the y- 
formalism (see Sect. [3}. Typically, the first CE is described by 
the y-scenario and the second by the a-formalism, if mass trans- 
fer is unstable. 

We applied the updated version of the population synthesis 
code SeBa to simulate the population of DWDs and SNIa pro- 
genitors. At present, close DWDs (of all WD types) are the clos- 
est related systems to the DD SNIa progenitors that are visible 
in bulk. The mass ratio distribution of the DWDs in model aa 
is inconsistent with the observations. Using model ya the sim- 
ulated population of DWDs compares well with observations, 
nevertheless, t his is wha t the y-formalism was designed to do. 

Recently, IWebbinkl d2008l) and IZorotovic etaD d2010h 
claimed that the predictive power of the ^-scenario is more 
restricted. They suggested that the a-scenario is valid when 
sources of an energy other than the binding energy of the en- 
velope is available, such as, the energy released by recombina- 
tion in the common envelope. This could explain the high value 
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Fig. 8. Simulated distribution of the population of merging double CO white dwarfs from a single burst of star formation as a 
function of delay time and total mass of the system. On the left model ya is used for the common envelope parametrisation, on the 
right model aa (see Sect. [3}. The intensity of the grey scale corresponds to the density of objects on a linear scale in units of number 
of systems per 1O 5 M . The black line corresponds to a combined mass of 1.4M . 
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Fig. 9. Simulated population of merging double CO white dwarfs from a single burst of star formation as a function of the masses 
of the two white dwarfs. M mass i ve is the mass of the most massive white dwarf, M non _ mass j ve corresponds to the least massive white 
dwarf. On the left model ya is used for the common envelope parametrisation, on the right model aa (see Sect. [3}. The intensity 
of the grey scale corresponds to the density of objects on a linear scale in units of number of systems per 10 s M . To increase the 
contrast, we placed an upper limit on the intensity, which effects only one bin for model ya and two bins for model aa. The black 
solid line corresponds to a combined mass of 1.4M . The dashed and dashed-dotted line correspond to a mass ratio q = m/ M of 
1 and 0.8, respectively. Th e mass ratio distribution, which is important e.g. in the scenarios proposed by Pak mor et alJ (j2010) and 
Ivan Kerkwiik et ail d2010h are very different for model ya and aa. 



of a found by iNelemans et al.l {2000) for the second CE, but 
certainly does not solv e the problem for the first CE for which 
INelemans etafld2000l) found a value of a < 0. 

The delay time distributions from our two models show the 
characteristic shape of a strong decay with time. This strong 
decay is expected when the delay time is dominated by the 
gravitational wave timescale (t gr °c a 4 ) and the distribution 
of orbital separations at DWD formation is si milar to th e ini- 
tial (ZAMS) distribution of N(a)da oc cC da (lAbtlll983h . The 
DTD from model ya fits t he shape of the observed DTD 
best. iMennekens et al.l {2010) also snowed a DTD using the y- 
scenario for the CE phase. They found that the DD DTD lies 
almost an order of magnitude lower in absolute rate than when 
using the a-scenario. However, they used the y-formalism for 



all CE phases. In our prescription (see Sect. [3]) the y-formalism 
is typically used in the first CE phase only. The reason for 
this is that in equal mass systems there is more angular mo- 
mentum compared to un equal mass systems with simila r orbits . 
I Mennekens et al.l d201C) and also Yungelson & Liviol {2000); 
kuiter et all d2009t) and lciaevs et alJ d201 lh showed DD DTDs 
using the a-scenario (however their CE-prescriptions may dif- 
fer slightly from Eqj2]>. Surprisingly, but as realised before, even 
though different groups used different binary evolution codes 
with different versions of the a-CE and CE efficiencies, the 
DTDs of the DD channel are very unifor m in that they show a 
strong decline with time (see for example Nele mans et al.ll20T2l 
for an overview). 
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Usually in synthetic DTD studies, the shape and normali- 
sation of the DTD are discussed separately. This might not be 
valid any more, as more and more observed rates are available 
and the conversion from observational units to synthetic units 
(e.g. the star formation history (SFH) and rate in per K-band lu- 
minosity instead of per M of created stars) is better understood. 
For example, the SFH is often convolved with the DTD to es- 
timate the SNIa rate in spiral galaxies like our Milky Way. The 
problem with this is that different assumptions for the Galactic 
SFH can significantly alter the theoretical Galactic SNIa rate. 
Since the SNIa rate follows the SFH with typical delay times 
of a few Gyr, the synthetic Galactic SNIa rate is very sensi- 
tive to the assumed SFH at recent times. When a constant SFH 
(of > 3M yr _1 ) is assumed, the SNIa rate is artificially en- 
hanced compared with detailed SFHs that show a peak in the 
star formation at a few Gyr and a de cline to lM Q yr _I at recent 
ti mes, see e .g. Nelemans et al. (2004). I n the observ ed SNIa rates 
of iMaoz & Badenesl (120101) and lMaoz et al.1 d2011l) the detailed 
SFH of every individual galaxy or galaxy subunit was taken into 
account to reconstruct the DTD. Therefore it is no longer neces- 
sary to convolve the theoretical SNIa rate from a burst of star for- 
mation with an approximate SFH. The theoretical calculations 
of the SNIa rate from a single starburst can directly be compared 
with observations. 

We found that the normalisation of the DTD of model aa 
and ya do not differ much, even though the CE evolution is very 
different. The time-integrated number of SNIa in model aa (3.3 • 
IO^Mq 1 ) is 70% larger as in model ya (2.0- IO^Mq 1 ). But most 
importantly, the simulated time-integrated numb ers do not match 
the observed number of 2.3±0.6- lO^Mg 1 by dMaoz etal.l201ll) 
by a factor of ~ 7 - 12. If our understanding of binary evolution 
and initial parameter distributions is correct, the standard DD 
channel is not a major contributor to the SNIa rate. 

For the SNIa model proposed by iPakmor et al ] (120101) . in 
which carbon burning is ignited in the merger process of two 
massive white dwarfs of nearly equal mass, we found an SNIa 
rate of 2.3 ■ lO^M' 1 for model ya and 1 .8 ■ 1O~ 5 M ' for model 
aa. IPakmor et al.l (1201 Oh founds that these systems resemble 
subluminous SNIa such as SN 1991bg. Assuming the fraction 
of 1991bg-l ike events to all SNIa events is 15+6% dLi et al.1 
200 ll 1201 ll), the o b served event rate is 3.5 ± 1.6 • lO^M" 1 . 
van Kerkwijk e t al.l (1201 Ol) proposed a model in which sub- 



Chandrasekhar WDs can explode as an SNIa. In this scenario 
two white dwarfs of nearly equal mass merge, though carbon ig- 
nition occurs only after the merger when the thick disk surround- 
ing the remnant is accreted onto it. The event rate is 2.3-lO~ 4 M 4 
for model ya and 1 .9 • lO^M" 1 for model aa. When only taking 
into account systems with a combined mass below 1 AM Q , the 
rates are 1.3 • lO^M ' 1 and 8.8 ■ 10~ 5 M~j , respectively. In the 
sc enario proposed by IPakmor et al.l (1201 Oh and in the scenario 
bv lvan Kerkwiik et al.l (l2010l) . systems are required to have high 
mass ratios. We showed that the mass ratio distribution of DWDs 
depends on the prescription for the CE. When the y-scenario is 
used, the average mass ratio of DWDs lies closer to one, which 
increases the SNIa rate in the above described scenario with re- 
spect to the a-scenario. T he rates of the channel proposed by 
Ivan Kerkwiik et ail d2010h for systems with sub-Chandrasekhar 
masses are on the same order of magnitude as the rates of the 
standard DWD channel. Therefore the combination of the two 
models is not sufficient to explain the observed rates. For our 
synthetic rates of the DD scenario to match the observed SNIa 
rates, within our current model of binary evolution, the param- 
eter space of the DD progenitors has to be increased severely, 



e.g. to include all CO-CO and CO-He mergers, which seems un- 
likely. 

Alternatively (and if contributions from channels other than 
the DD are minor), our model underpredicts the fraction of stan- 
dard DD SNIa progenitors in the entire DWD population. Our 
model of the visible population of DWDs predicts 0.9-2.9% of 
the visible DWDs (depending on th e model) to b e SNIa progen- 
itors. To match the observed rate of IMaoz et al.l d201 ll) . 10-30% 
(excluding any errors on the observed and synthetic values) of 
the observed DWDs should lie in the SNIa progenitor region 
(upper left corner of Fig. [3). With 46 observed DWDs so far, 4- 
15 SNIa progenitors are expected without taking non-uniform 
selection effects into account. So far, only two systems have 
been found that possibly are SNIa progenitors, which makes it 
improbable, but not impossible, that our model underpredicts 
the number of DD SNIa progenitors. When the population of 
observed DWDs is increased, the fraction of SNIa progenitors 
amongst DWDs will give more insight into the validity of our 
knowledge of binary evolution of massive DWDs. 

Concluding, although the shape of the DD DTD fits the ob- 
served DTD beautifully, the normalisation does not. An impor- 
tant point is that we did not optimise our model to fit the ob- 
served DTD in shape or number. We showed that the normal- 
isation can be influenced by the metallicity; ~ 30 - 60% de- 
pending on the model for Z = 0.001 with respect to Z = 0.02. 
Furthermore, the normalisation depends on the initial mass func- 
tion, the percentage of single stars, and the initial distribu- 
tion of mass ra t ios an d orbital periods. In this paper and in 
iNelemans et al.l (1201 2|) we assume d the percentage of single 
star s to be 50%. Results from e.g. iKouwenhoven et all (120071) 
and lRaghavan et aLl d2010l) showed that the binary fraction might 
be as high as 70% or more for A- and B-type stars, potentially 
raising the synthetic SNIa rate by a factor < 2. Preliminary 
results show that the initial distribution of mass ratio and or- 
bital separation affects the slope of the DTD, still the strong 
decline with time remains. Moreover, the integrated rates are 
not affected by factors sufficient to match the observed rate. 
Additional research is needed to study if the normalisation can 
be raised sufficiently to match the observed rate. If not, the 
main contribution to the SNIa rate comes from other channels, 
such as the SD scenario (e.g. supersoft sour ces), double detonat - 
ing sub-Chandrasekhar accretors (see e. g. iRromer et al.l 20101). 
or Ko zai oscillations in triple systems (Shappee & Thompson 
1201 2t Hamers et al. in prep.). 
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Appendix A: Most important changes to the 
population synthesis code SeBa 

A.1. Treatment of wind mass-loss 

Each star may lose mass in the form of a stellar wind. In a bi- 
nary system the stellar wind matter from a binary component 
can be accreted by the companion star or lost from the system 
(see Appendix IA.2I ). This influences the binary parameters via 
the loss of mass and angular momentum from the system. We 
assume that the matter that is lost from the system carries a spe- 
cific angular momentum equal to that of the star from which 
it originates. Furthermore, we assume that wind accretion onto 



the bi nary companion is Bon di-Hoyle accretio n ( Bondi & Hovld 
1 1944 . as re-formulated by iLivio & Warned (11984 . The wind 
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mass loss prescriptions for different types of stars used in SeBa 
are updated e.g. to include metallicity dependency where possi- 
ble. The prescr iptions correspond t o some degree to the recom- 
mendations by iHurlev et al.l d2000l) . If multiple mass loss pre- 
dictions are applicable to a star, we take the one that predicts the 
maximum mass loss rate. 

- For all types of luminous stars (L > 4OOOL ) from 
the main sequence (MS) to the asymptotic giant branch 
(AGB) we apply the empiric al mass loss rate by 
iNieuwenhuiizen & de Jageri (Qj?90) given by 

M = 9.631 ■ 1(T 15 R (YU L 124 M 016 M Q yr -1 , (A.l) 

where M is the mass accretion rate, R the stellar radius in 
R Q , L the luminosity in L and M the stellar mass in M Q . We 
assum e that the formalism of INieuwenhuiizen & de Jageri 
(1990) is dependen t on the initial meta llicity as M(z) = 
(z/z ) 1/2 M(z ) (see iKudritzki rtlill987h . 

- For a mas s ive MS star we give preference to the rates of 
IVink et al.1 d2000l 1200 lb. Where they do not apply, the rates 
of Nieuwenhuiizen & de Jager ( 1990) are used. Massive MS 
suffer from strong win ds driven by radiation p ressure in lines 
and in the continuum. IVink et al.l A 20032001) take into ac- 
count multiple scattering effects of photons. They find good 
agreement between observations and theoretical mass-loss 
rates. 

- For stars in giant phase s we adopt the empirical relation 
found bv iReimersI (119751) . 

M = 4- 10~ 13 — M yr 1 . (A.2) 
M 

We assume a numeric al p refactor of n _ 0.5, see 
lMaeder&Mevnetl(ll989h and lHurlev eUdlfcoOOl) . 

- AGB stars can experience severe mass-loss caused by radi- 
ation pressure on dust that condensates in the upper atmo- 
sphere of the stars. Empirically, the mass-loss rate has been 
coupl ed to the period of large- amplitude radial pulsations 
Ppuis (IVassiliadis & Woodlll993l) : 

log P puls (days) = -2.07 + 1 .94 ■ \og(R) - 0.9 • log(M). (A.3) 

We apply mas s-loss to the enve l ope ac cording to the pre- 
scription of IVassiliadis & Wood! dl993l) . During the super- 
wind phase the radiation pressure driven wind is modelled 
by 

M = , (A.4) 

^exp 

where c represents the speed of light and v exp the stellar wind 
expansion velocity. The latter is given by: 

v exp (km s- 1 ) = -13.5 + 0.056P pu i s (days). (A.5) 

Furthermore, v exp is constrained to the range 3.0-15.0 km 
s" 1 . 

Before the superwind phase, the mass loss rate increases ex- 
ponentially with Pp U i s as 

logM(M yr 1 ) = (A.6) 

-11.4 + 0.0123 -Ppuis ifM<2.5, 
-11.4 + 0.0125- (Pp U i s - 100 ■ (M - 2.5)) ifM > 2.5. 

The mass loss rate of Vassiliadis & Wood ( 1993) is given by 
the minimum of Eq. IA.4l and lA.6l 



- Luminous blue variables (LBVs) are extremely massive 
and luminous stars near the Humphreys-Davidson limit 
dHumphrevs & Davidson! 1 19941) with enormous mass-loss 
rates. We use the L BV mass loss prescri ption and implemen- 
tation suggested bv lHurlev et al.l (Hp00): 

M = 0.1 X ( 10~ 5 RL l/2 - 1.0? ( ^— - l.o) M yr 1 , 

v ' \6.0- 10 5 / 3 

(A.7) 

if L > 6.0 ■ 1O 5 L and 10^L 1/2 > 1.0. 

- Wolf-Rayet stars are stars in a stage of evolution following 
the LBV phase where wea k or no hydrogen li nes are ob- 
served in their spectra. Like lHurlev et al.l d2000h we include 
a Wolf-Rayet-like mass-loss for stars with a small hydrogen- 
envelope mass (ji < 1.0 from their Eq. 97). The prescrip- 
tion itself, however, is different. W e model it according to 
iNelemans & van den Heuvelld200l : 

M = 1.38 • 10~ 8 M 2 - 87 M Q yr 1 . (A.8) 

This is a fit to observed mass-loss rates from 

Nugis & L amersl ([2000). We multiply with a factor 
(1 - jj) to smoothly switch on mass loss. 

- In addition to the evolution of ordinary hydrogen rich stars, 
the evolution of helium burning stars with hydrogen poor 
envelopes is simulated as well. For helium main-sequence 
stars with a mass M > 2.5M we assume the same relation 
as for Wolf-Rayet-like stars. For helium giants, either on the 
equivalent of the Hertzsprung or giant bran ch, we describe 
mass lo ss in a very general way similar to INelemans et al.l 
(2001b). We presume 30% of the mass of the envelope M env 
will be lost during the naked helium giant phase with a rate 
that increases in time according to 

AM wmd = 0.3M env [(— ) 6 ' 8 - (lf% (A.9) 
ff ff 

where AM w ; n d is the amount of mass lost in the wind in M in 
a timestep At, ff is the duration of the helium giant phase and 
t the time since the beginning of the phase. 

Special attention has been given to prevent large wind mass 
losses in single timesteps because the mass loss prescriptions 
are very dependent on the stellar parameters of that timestep. 
For this reason we implemented an adaptive timestep in situa- 
tions where strong winds are expected, e.g. at the tip of the giant 
branch. This procedure is accurate to differences in stellar mass 
of less than 4% for masses below 12M . 

A.2. Accretion onto stellar objects 

Roche lobe overflow mass transfer and subsequent accretion can 
substantially alter the stars and the binary orbit. Mass accretion 
can affect the structure of the receiver star and its subsequent 
evolution. When more mass is transferred than the accretor can 
accrete, we assume that the non-accreted matter leaves the sys- 
tem with an angular moment um of 2.5 times the specific angu- 
lar momentum of the b inary dPortegies Zwart & Verbunli ll996: 
INelemans et al.ll2~001bl) . For compact accretors we assume the 
matter leaves the system with the specific angular momentum of 
the compact remnant. In this section we discuss the limiting ac- 
cretion rate, the response of the accretor to regain equilibrium, 
and the subsequent evolution of the new object for different types 
of accretors. 



13 



S. Toonen et al.: Supernova Type la progenitors from merging double white dwarfs: 



A.2.1 . Accretion onto ordinary stars 

For ordinary stars from MS to AGB stars, we distinguish be- 
tween two types of accretion; accretion from a hydrogen-rich or 
a helium-rich envelope. Hydrogen-rich accretion can occur for 
example when a donor star ascends the giant branch and fills its 
Roche lobe. After it loses its hydrogen envelope, it can become a 
helium-burning core. When this helium star ascends the helium 
equivalent of the giant branch, a fraction of the helium-rich en- 
velope can be transferred onto the accretor. We name this type 
of accretion 'helium accretion'. We assume that the accreted he- 
lium settles and sinks to the core instantaneously. The helium 
accretion rate is limited by the Eddington limit. Hydrogen is ac- 
creted onto the envelope of the receiver star. The accretion rate 
is bounded by the star's thermal timescale times a factor that is 
dependent of the ratio of Roche lob e radius of the receivers to 
its effe ctive radius, as described by Porteg ies Zwart & Verbunj 
(I1996I) . The formalism is proposed by Pols & Marinus (1994), 
which is based on Kippenhahn & Meyer-Hofmeistej] (1 19771) ; 
iNeo etalJ (1 19771) and iPacket & De Grevd (1 19791) . If the mass 
transfer rate is higher than the maximum mass accretion rate, 
the excess material is assumed to leave the binary system. 

Because of the accretion, the star falls temporarily out of 
thermal equilibrium. While regaining equilibrium, the gas en- 
velope surrounding the core puffs outward. Because we do 
not solve the equations of stellar structure and the stellar evo- 
lution tracks describe single stars in equilibrium, we add a 
procedure to account for a temp oral increase in radius as in 
Porte gies Zwart & Verbuntl (1996). This is important for exam- 
ple to determine if an accretor star fills its Roche lobe. It also af- 
fects the magnetic braking process and the Darwin-Riemann in- 
stability through the increased stellar angular momentum. Note 
that the mass transfer rate is not dependent on the stellar radius 
in our simulations, so that the binary evolution is not critically 
dependent on out-of-equilibrium parameter values. 

Accretion can also affect the structure of the receiver star and 
its subsequent evolution. It is modelled by changing the stellar 
track and moving along the track. The former is described by the 
track mass, which is equivalent to the zero-age main-sequence 
mass that the star would have had without interaction. The latter 
is described by the relative age t K \ of the star. We distinguish two 
cases: 

- Rejuvenation of an MS star 

Accretion onto an MS star rejuvenates the star. The star 
evolves similarly to a younger star of its new mass and its MS 
lifetime can be extended. It would show up in a Hertzsprung- 
Russell diagram as a blue straggler. For hydrogen accretion 
the track mass is always updated and the renewed relative 
age of the star 



l rel 



?rel 



t ms M' 



(A. 10) 



where primes denote quantities after a small amount of mass 
accretion, f ms the main-sequence lifetime, and M the mass of 
the star. 

For helium accretion we assume the mass accretes to the core 
instantaneously and the track mass is increased accordingly. 
These stars appear older than for hydrogen accretion because 
more hydrogen has been burned previously. The rejuvenation 
process is described by 



where we assume 10% of the mass of the star will be burned 
during the MS phase. 
- Rejuvenation of a giant 

During the giant phases the envelope is discoupled from the 
core in terms of stellar structure. The evolution of the star 
will therefore not be influenced directly by small amounts 
of hydrogen accretion to the envelope. The track is only up- 
dated when the new mass is larger than the track mass to 
account for severe hydrogen accretion. The mass before ac- 
cretion can be much lower than the track mass because of 
wind mass loss, which can be strong for giants. For helium 
accretion to the core, the track is always updated. An excep- 
tion to this is the early AGB where the helium core does not 
grow. In this stage there is a one-to-one relation between the 
helium core mass and the track mass (Eq.66 in lHurlev et al.l 
120001) . 

When a giant accretor star moves to a new evolutionary 
track, we need to determine the location of the star along 
this track. In a more physical picture this means determin- 
ing the relative age of the star t re \. For a giant its evolution 
is mainly determined by its core. Therefore for a given evo- 
lutionary track and core mass, the relative age is effectively 
constrained. For both types of accretion, we insist that the 
star stays in its same evolutionary state after its mass in- 
crease. When no solution can be found for ? le i, the relative 
age is set to the beginning or end of the current evolutionary 
state and the track mass is varied to find a fitting track that 
ensures mass conservation. 

A.2.2. Accretion onto helium-burning cores 

For accretion onto helium-burning stars that have lost their hy- 
drogen envelopes, accretion is limited by the Eddington limit. 
Helium accretion onto a helium main-sequence star is similar as 
hydrogen accretion onto normal main-sequences stars. We as- 
sume that the star evolves similarly to a younger star of its new 
mass according to Eq. IA. 101 where f ms should be replaced by 
the helium main-sequence life time. We assume that for helium 
giants the envelope is discoupled from the core in terms of stel- 
lar structure, as with hydrogen-rich giants. Therefore we assume 
that the evolution of the giant is not affected and only update the 
track when the new mass is larger than the track mass. 

The effect of hydrogen accretion onto helium stars is more 
complicated. If the hydrogen layer is sufficiently thick, the layer 
can ignite. This can significantly increase the radius of the star 
and essentially turn it into a born-again star on the horizontal 
or asymptotic giant branch. We studied the effect of hydrogen 
accretion to helium stars with stellar models simulated by the 
stellar evolution code STARS. This code models stellar struc- 
ture and evolution in detail b y solving t he stel lar structure equa- 



, 4s M 6Mt' ms 

rel rel f ms M' 0AM'' 



(All) 



tions. The code is based on Eggleton ( 197 ~D) and includes up 
dated input physics as described in iHuetalJfcoich . The models 
do not include atomic diffusion. For mass accretion rates at ten 
percent of the Eddington rate of the accretor for ~ 10 4 - 10 5 yr, 
the accreted hydrogen layer ignites. Helium stars that are more 
massive than ~ 0.55M Q resemble horizontal branch stars after 
accretion, but most of the luminosity still comes from helium 
burning. For lower mass helium stars this is not the case, be- 
cause the corresponding horizontal branch stars (< 3.5M©) can 
have ignited helium in a degenerate core, which strongly affects 
the characteristics of the star. For both mass ranges, the accre- 
tor expands by a factor ~ 10-100 compared to the original 
helium star. Because hydrogen accretion to helium stars is not 
very likely, we model this very simply. When more than 5% of 
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the total mass is accreted, the radius of the star is increased by a 
factor 50. With few exceptions, this leads to a merger of the two 
components. 

The effect of hydrogen accretion to helium giants is not 
known very well and additional research is necessary. For now, 
because it is very unlikely to happen, we treat it in the same way 
as helium accretion onto the envelope of the giant. 

A.2.3. Accretion onto remnants 

White dwarf, neutron star and black hole accretors can accrete 
with a maximum rate of the Eddington limit. If more mass is 
transferred, the surplus material leaves the system with the spe- 
cific angular momentum of the compact remnant. For neutron 
stars and black holes we assume that the transferred mass is tem- 
porarily stored in a disk. From this disk, mass will flow onto the 
surface of the remnant with ten percent of the Eddington limit. 
We assume that a neutron star collapses onto a black hole when 
its mass exceeds 1.5M Q . 

For white dwarfs, the accretion process is more complicated 
because of possible thermonuclear runaways in the accreted ma- 
terial on the surface of the white dwarf. In SeBa there are sev- 
eral options to model the effectiveness of the white dwarf to 
retain the transferred material. For hy drogen accre t ion w e can 
choose between the efficiencies of Hachisu et al ] (120081) and 
iPrialnk &JCovet^([l995 [). For helium retention, th e option is be- 
tween Kato & Hachisu (fl999l) (with updates from lHachisu et all 



119991) and llben & TutukovT (ri996) 



A3. Stability of mass transfer 

A semi-detached system can become unstable in two ways. In 
a mass transfer instability, the Roche-lobe-filling star expands 
faster than the Roche lobe itself on the relevant ti mescale. In th e 
other case tidal interactions lead to an instability (iDarwi n 1879). 



A.3.1. Tidal instability 

A tidal instability can take place in systems of extreme mass ra- 
tios. When there is insufficient orbital angular momentum 7b that 
can be transferred onto the mass-losing star, the star cannot stay 
in synchronous rotation. Tidal forces will cause the companion 
to spiral into the envelope of the donor star. The tidal instability 
occurs when the angular momentum of the star 7+ > I 7b, where 



J h = Mm 



Ga 
M + m 



(A. 12) 



where 7b is the orbital angular momentum of the circularised 
binary, a is the orbital separation, M the mass of the donor star 
and m the mass of the accretor star. The angular momentum 7* 
of a star with radius R is given by 



7. 



k 2 MR 2 a>, 



(A.13) 



where k 2 is the gyration radius described by iNelemans et al.l 
d2001bl) and a> is the angular velocity of the donor star, which 
is assumed to be synchronised with the orbit. It is given by 
a) = 2n/Pb, where P\, is the orbital period. We model the in- 
spiral according to the standard a-CE (see Sect. 0. Owing to 
the expulsion of the envelope, the binary may evolve to a more 
stable configuration or merge. If the mass-losing star is a main- 
sequence star, we assume that the instability always leads to a 
merger. 



A.3.2. Mass transfer instability 

The stability of mass transfer from Roche lobe overflow and its 
consequences on the binary depend on the response of the ra- 
dius and t he Roche lobe of the donor star to the impo sed mass 
loss (e.g. Web binklll985UHieUming & Webbinkfll 9871 (hereafter 
HW87); iPols & Marinuslll994t ISoberman et al.lfl997l) . We dis- 
tinguish four modes of mass transfer; on the dynamical, thermal, 
nuclear timescale of the donor or on the angular-momentum-loss 
timescale. The response of the accretor star to the mass that is 
transferred onto it and the effect of this on the orbit is described 
in Appendix lA.2l The response of the donor star to mass loss is 
to readjust its structure to recover hydrostatic and thermal equi- 
librium. The dynamical timescale to recover hydrostatic equilib- 
rium is short compared to the thermal timescale. For mass trans- 
fer to be dynamically stable, the dynamical timescale of the star 
is important. The change in radius due to adiabatic adjustment of 
hydrostatic equilibrium is expressed as a logarithmic derivative 
of the radius with respect of mass, 



fad 



d\nR 
d\nM 



(A. 14) 



where M and R are the mass and radius of the donor star. The 
assumed values of f a d are shown in Table IA.il 

The response of the Roche lobe R^ of the donor star is ex- 
pressed as the logarithmic derivative of the Roche lobe radius 
with respect to mass: 



d\n R L 
dlnW 



(A.15) 



The value of £l is calculated numerically by transferring a test 
mass of 10~ 5 M . Because £l = £iXM\, Mi,d), fx is dependent 
on the mass accretion efficiency of the secondary, and therefore 
on the mass accretion rate of the test mass. For instance, for high 
mass ratios q » 1 the loss of some mass and corresponding 
angular momentum can have a stabilising effect on the mass- 
transferring binary. To determine the dynamical stability of mass 
transfer, we assume that the mass transfer rate of the test mass is 
on the thermal timescale of the donor star: 



GM 2 
RL 



(A. 16) 



1. When £l > fad, mass transfer is dynamically unstable. We 
model this as a CE phase, as described in Sect. [3] 

When £l < fad mass transfer is dynamically stable. The 
donor star is able to regain hydrostatic equilibrium and shrinks 
within its Roche lobe on a dynamical timescale. To determine if 
the donor star is also able to regain thermal equilibrium during 
mass transfer, we calculate the change in the radius of the star as 
it adjusts to the new thermal equilibrium: 



feq — 



<ilnR 
d\nU 



(A. 18) 



2 When the timescale of the CE phase becomes relevant, we as- 
sume that it proceeds on a time scale r given by the geometric mean 
of the thermal T t h and dynam ical Td timescales of the donor (see 
Paczyinski & Sienkiewicz 1972): 



(A. 17) 
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The assumed values for £ eq are described in Table IA. 11 

To calculate the response of the Roche lobe ^L,eq, we assume 
that the mass transfer rate of the test mass is on the nuclear evo- 
lution timescale: 

T mc = dt—-, (A. 19) 

where R represents the equilibrium radius of the star according to 
the single-star tracks. dR eq is the change in R in a short timestep 
dt without binary interactions. 

2. When ^L,eq < nw((eqi <Tad) mass transfer is driven by the ex- 
pansion of the stellar radius due to its internal evolution. 

3. When ^L,eq > '«'«(£eq, £adX the mass transfer is thermally 
unstable and proceeds on the thermal time scale of the donor. 

4. The previous modes of mass transfer are caused by an ex- 
panding donor star. The final mode is caused by shrinking of 
the orbit caused by angular momentum loss. We assume that 
this mode takes place when the corresponding timescale Tj is 
shorter than the timescales at which the other three modes of 
mass transfer take place. Angular momen tum loss can hap- 
pen due to gravitational wave radiation ■L- (iKraft et all 19621) 

I III I II I 

and magnetic braking J m \, (Schatzman 1962; Huang 1966; 
ISkumanichl Il972t IVerbunt & Zwaanl 1 19811) Mass transfer 
proceeds on the time scale on which these processes occur: 



Jb 



J QT + J m \y 



(A.20) 



where 7b is the angular momentum of the circularized binary 
given by Eq. IA.12l Next we discuss the assumptions and im- 
plications of igr and / m b. 

Gravitational wave radiation most strongly influences close 
binaries since it is a strong function of orbital separation. 
The change in orbital sep aration d averaged over a full orbit 
is given by dPeterslll964l) 



64 G 3 Mm(M + m) 
~~5 c 5 a 3 (l - e 2 ) 7 - 2 



where j gl -/Jb = a gr /(2a). 

Magnetic braking mostly affects low-mass stars within the 
mass range of 0.6 < M/M Q < 1.5. These stars suffer from 
winds that are magnetically coupled to the star. Although the 
mass loss in this process is negligi ble, the associated a ngular 
momentum loss can be severe (Rappapo rt et al.lll983b : 



J mb 



3.8 ■ KT 30 /^ P {M + m)RP 



0} 



ma- 



(A.22) 



where /3 is a parameter that represents the dependence of the 
braking on the radius of the donor star. We take f3 — 2.5. 



Table A.l. Values of the adiabatic £ a( j and thermal £ eq response 
of the radius to mass loss for different types of stars. 



k 


Evolutionary state 






0,1 


Main-sequence: 








M a < 0.4 


-1/3 







0.4 <M a < 1.5 


2 


0.9 




M a > 1.5 


4 


0.55 


2 


Hertzsprung gap: 








M a < 0.4 


4 







M a > 0.4 


4 


-2 


3 


First giant branch: 








■ shallow convective layer 


4 







• deep convective layer 


HW87 





4 


Horizontal branch: 








M a < M Hcf 


4 


4 




A? Hc f < M a < M FGB : 








■ decent along GB 


as k = 3 







■ blue phase 


4 


4 




M > M FGB : 








■ blue phase 


4 


-2 




• ascent to AGB 


HW87 





5,6 


Asymptotic giant branch 


HW87 





7 


Helium main- sequence: 








M < 0.2 


15 


-0.19 




M > 0.2 


15 


1 


8,9 


Helium giant: 








M c < 0.4 


HW87 


-1/3 




M c > 0.4 


HW87 


-2 


10,11,12 


White dwarf 


-1/3 


-1/3 



Notes. The types of stars correspond to the definition by Hurle y et al.1 
(2000) expressed by their integer k. The stellar tracks are distinguished 
by the mass M a , which is equivalent to the ZAMS mass the star would 
have had without interac tion. Mn c f and Mfgb are defined by Eq. 2 and 3 
from lHurlev et alJd2000l) . M Hc f represents the maximum initial mass for 
which helium ignites degenerately in a helium flash, which is ~ 2M for 
solar metallicities. M FG b is the maximum initial mass for which helium 
ignites on the first giant branch, which is ~ 13M Q for solar metallicities. 
M is t he total mass of the star and M c the mass of the core. HW87 rep- 
resents Hiellming & Webbink| d!987l) . who calculated gad f° r condensed 
polytropes, consisting of a compact core surrounded by an envelope 
with polytropic index n = 3/2. For stars on the first giant branch there 
are two prescriptions of £,d. If the convective zone in the upper layers of 
the envelope is shallow (fits from Yungelson, private communication), 
we assumed the envelope responds to mass loss in a similar manner as 
radiative envelopes. 
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